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I. INTRODUCTION 


The origin of high energy cosmic rays (CRs) has been a long-standing problem in astrophysics. Especially, CRs with 
energy above > 10^® eV are considered to come from extra-Galactic sources such as active galactic nuclei (AGNs) 
and gamma-ray bursts (GRBs). In these sources we expect the production of high energy neutrinos (> 0.1 TeV) 
through interactions of accelerated protons with the ambient photons (py interactions) or gas {pp or pn interactions) 
GHi. Detection of these neutrinos can provide us new information about high energy cosmic ray sources as well as 
the acceleration processes. 

In cosmic ray accelerators, high energy neutrinos are mainly produced from the decay of charged pions: tt®' —>■ 
^ + Vfj, + Ve + and 7r“ —>■ p~ + D ^ e~ + Therefore, the flavor ratios of these neutrinos 

are expected to be 


$0 : = 1 : 2 : 0 , 

t/g Un. ' 


( 1 ) 


at the sources, where denotes the flux of Vi and Vi {i = e, p or r). The observed flavor ratios become : 

=1:1:1 after the neutrino oscillations during the propagation to the Earth However, this argument may be 
too naive because we should take into account the finiteness of the decay timescale of pions and muons p^. For 
example, if the cooling timescale or acceleration timescale of a pion or a muon is shorter than the decay 

timescale, the spectral shape of neutrinos produced from the decay of those particles would be significantly modified. 
Especially, because the decay times are different between pions and muons, the energy dependence of neutrino fluxes 
would be different from flavor to flavor. The observed neutrino flavor ratio may be also modified by neutron decay 
and new physics such as neutrino decay Gl, [l^ , sterile neutrinos G^ i pseudo-Dirac neutrinos SGI , Lorentz 
or CPT violation G3> quantum-gravity G3 secret interactions of neutrinos [^ . 

Recently the high energy neutrino detector, IceCube, has discovered 30 TeV — 2 PeV neutrinos 2l| which are 


confirmed to be non-atmospheric at the level of 5.Ter. The IceGube team has also analyzed the flavor composition of 
astrophysical neutrinos in the energy range of 35 TeV — 1.9 PeV and demonstrate consistency with = 

1:1:1 (although the best fit composition is 0 : 0.2 : 0.8). In the near future, the next generation IceGube-Gen2 and 
KM3Net experiments will enable the precise study of the energy spectrum of high energy neutrinos and their flavor 
composition. 

There is no study about how the flavor ratio of observed neutrinos as well as its energy dependence are modified 
by the acceleration of pions and muons, although several authors have investigated the neutrino energy spectrum 
under the secondary-acceleration iSil. In addition, their approaches are based on one-zone [1, [l^ or two-zone Ga 
approximations, that is, they do not consider the spatial distribution of secondary particles (pions and muons) and 
their transport across the shock. 

In this work, we investigate the acceleration of pions and muons produced by protons by solving their convection- 
diffusion equations around the shock front taking into account their decay into other particles (i.e., pions into muons 
and muon neutrinos, and muons into muon neutrinos and electron neutrinos). The shock acceleration of secondary 
particles has been discussed in the context of the positron excess [22l - [^ observed by PAMELA/Fermi LAT/AMS-02 
pbl - l^ (see also [2§-[3l|). We develop their formalism by including the decay of secondary particles during their 


acceleration, and evaluate the energy spectrum of neutrinos produced from those particles. 

This paper is organized as follows. In Section 2 we formulate a general model for the shock acceleration of pions and 
muons, which are produced from shock-accelerated protons via photomeson interactions, and give versatile expressions 
of the neutrino spectra for later application. In Section 3, we show that the acceleration of pions/muons would be 
possible in low-power GRBs occurring inside their progenitors, and apply our model to that case to compute the 
neutrino spectra and flavor ratios. Our results and discussions are summarized in Section 4. In Appendix A, we 
summarize general solutions of the convection-diffusion equations for pions (decaying secondary particles) and in 
Appendix B for muons (decaying tertiary particles). 


II. MODEL 


In this section, we describe the shock acceleration of secondary pions and muons that are generated from the 
protons accelerated at the shock, and investigate the energy spectrum of neutrinos produced by those pions and 
muons without specifying a particular source. Hereafter we neglect the energy loss of particles due to synchrotron 
emission and inverse Gompton scattering for simplicity (see discussions in Sec. 4). 

In the shock rest frame, the transport and shock acceleration of particles decaying into other kinds of particles with 
a timescale can be described by the convection-diffusion equation, 



A 

dx 


nf 


P du dfi 
3 dx dp 


h 

Ti 


Q^{x,p), 


( 2 ) 





3 


where fi{x,p) {i = 7r,/r) is the equilibrium distribution function of accelerated particles per unit spatial volume and 
per unit volume in momentum space, D(p) is the diffusion coefficient, and u is the the velocity of the background 
fluid. We assume that the shock is non-relativistic and that the distribution functions are stationary and isotropic 
except for the shock front for simplicity. When the shock is relativistic, we should solve the relativistic version of 
the convection-diffusion equation taking into account the anisotropy of the particle momentum distribution. The 
anisotropy is the order of P-, and therefore in the mildly-relativistic shock, it would be less than a factor of two. 

The shock front is set at a; = 0, and the upstream (downstream) region corresponds to x < 0 (x > 0). If we ignore 
the third term of the right-hand side, which describes the loss of particles due to their decay, this is a well-known 
equation for the usual diffusive shock acceleration of particles . The decay timescales of a pion and a muon 

are the functions of their energy. 


'^TT — ^7r,0 " 




~ 1.9 X 10 ^ S £,r,l00TeV, 


= 'T/i.O- 


~ 2.1 S £^,100TeV, 


( 3 ) 

( 4 ) 


where Si = 100 TeVeipooTeV is the energy of a particle at the shock rest frame and, — 2.6 x 10“® s and 
— 2.2 X 10“® s are the decay timescales of a pion and a muon at their rest frames, respectively. 

The fourth term of the right-hand side of Eq.Q, Qi(x,p), is the distribution function of i particles injected per unit 
time. Here we consider that charged pions are produced in py interactions, pj —>■ A+ —>■ 7r+n. In this case, Qjr(x,p) 
should be given from the distribution function of primary protons. On the other hand, in the case of muons, since 
they are produced by the decay of pions, Q^{x,p) should be given from the distribution of pions (see below). 

Hereafter, we assume the Bohm-type diffusion, 


D{p) 


rjc^p 

1^" 


( 5 ) 


where e is the c harg e of a particle, B is the magnetic field strength and 77 is the gyrofactor, which is equal to unity in 
the Bohm limit [^. The fluid velocity is given by 


u{x) 


it_ {x < 0), 
n_|_ {x > 0), 


where U- and are constants and the compression ratio is cr = m_/u+ > 1. 

One should solve Eq.@ taking the following boundary conditions into account: 


(i) lim fi = lim /j, 

X—^ —0 X —>-+0 

(ii) lim fi = 0, lim /* < 00, 

X—>■ —00 X^ + CXD 


(iii) 


nf 


x=—0 


. x=+0 


1, . dfi 

= 3 («+ - n-)v 


x—0 


( 6 ) 

( 7 ) 

( 8 ) 

(9) 


where (iii) comes from the integration of Eq. m across the shock front. This condition yields the differential equation 
for the distribution function at the shock front fip{p) = fi{x = 0,p) with respect of p (see Appendix A). 

The detail of the general solution of Eq. ([^l) is presented in the Appendix A. In the following subsections we briefly 
describe the properties of the derived distribution functions of pions and muons. 


A. pion acceleration 

In this subsection we show the pion distribution function evaluated from Eq.([2]). Pions are produced through the 
interactions between the protons accelerated in the shock and the ambient photons. The distribution of protons is 
given by 

=Its (">S: (“) 

where fp o(p) is the proton distribution function at the shock front, which is proportional to ^ P~"^■ This expression 
m is a well-known solution of the convection-diffusion equation ([5]) . 
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For simplicity we assume that the ambient photon field is uniform. Since pions are produced from shock-accelerated 
protons, the production spectrum at the source of pions is proportional to that of primary protons, which can be 
described as 


n r-r- - /‘57r,o(p)exp[a;u_/£>(pp)] (x < 0), , . 

(ill 

where Pp is the momentum of a primary proton generating a secondary pion with a momentum p, and these momenta 
are approximately related in a linear way: 


P ~ CttPp) (12) 

where ~ 0.2 is the ratio of the energy of a pion to that of a primary proton [I|. 

We can solve Eq. @ for pions by using Eqs. (Tni) . (IMI) . and (IA3I) with Qt^{x,p) in Eq. (ITT]) , and obtain the pion 
distribution functions in the upstream /,r,-(x,p) and downstream /,r,-i-(x,p) as 


U,- = 


f-rrfi — 


DQ 


TT.O 


exp 


^ 7 r,+ (/tTjO Q 7 r, 0 ^ 7 r)^^P I 


d/t^ + {in -iiW- 

\Ju\+ ^DjTn - U+ 


v'lT7lD7i“+»- \ 


2D 


X I Qn,0■ 


, ( 13 ) 

(14) 


The pion distribution functions at the shock front = fn{x = 0,p) can be evaluated from Eqs. (IA5I) . (IA6I) and 
(|A7[) as 


Jo P \pj ut 


(15) 


where and are numerical factors, being independent of p (since both D and are proportional to p): 


A-n — — 
2 


Bn = 


'1-b 


AD 

TnvJ_ 


+ 1 + < / ^ + 


1 4£» 


9 ' 2 

(7 TtzU- (7 


2(7 


ADlTnV?_ - (1 - 2in) AD/Tnu\ + 1 


(16) 

(17) 


where a = it_/u+ is the compression ratio. We can see that the distribution function of pions at the shock front, 
Eq. (IT5|) . becomes harder than their production spectrum Qnfi by p^ [oc D{p)\, and this is similar to Eq.(6) of [^ . 
where the acceleration of secondary positrons produced in the supernova remnant shock is discussed. The difference 
is that we take into account the decay of particles, the third term of the right-hand side of Eq. m, in the convection- 
diffusion equation of secondary particles, which is reflected as numerical factors A^ and B^- 

To see the effects of the transport and acceleration of pions in the downstream region from Eq. (HI, we divide 
fn,+{x,p) into two components: fn,a.cc{x,p) and fn,nonacc{x,p). The former component fn,acc represents the pions that 
are reaccelerated at the shock, being proportional to D{p)Qn,o(j>)- On the other hand, the latter component fn,nonacc 
represents the pions that are produced from the protons and advected in the downstream region, being proportional 
to Qn,oip)Tn-. 


/7r,acc(x,p) = fn,o{p)exp 




+ ADfTn — U+ 



(18) 


/"TT.nonacc(x, p) — Q7r,o{p')'^7r 


1 — exp — 


\Ju\+ ADjxn - u+ 


2D 


(19) 


Hereafter, since the number of upstream pions, fn,-{x,p), is subdominant compared to that of downstream 

pions, fn,+ {x,p), we only discuss the contribution of the downstream pions to the neutrino spectra. 
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In the limit of D/tt^v?_ —>■ 0 (i.e., the lifetime of a pion is much longer than the acceleration timescale, tacc = D/v?_), 

1, i?_ Ri +(7. Therefore, the distribution 


Ea. dTSl) has the same form as Eq.(6) of [22|. In this limit, we have 
function at the shock front in Eq- m can be approximated as 


/tt.O — 


7 


7 - a + 1 VC, 


1 


D{p)Q^,oip) 


( 20 ) 


where a is the power-law index of the production spectrum, Qt^^Ip) oc p~°‘, and a r; 4 in the strong shock limit 
with an adiabatic index 5/3. We can see that, since we assume D{p) oc p, the resulting spectrum is proportional to 
being harder than the production spectrum. This can be interpreted as the result of the secondary-acceleration: 
since pions produced from shock-accelerated protons can cross the shock front before their decay, they would gain the 
energy and their spectrum would become harder. We should also note that is proportional to taccQ^.o- In this 
limit of DItt^u^ —>■ 0, the pion distribution function in the downstream region (1141) can be approximated as 


/77.+ ^ /7r,o(p)exp 


— ) 

U+T^ J 


^7r,o(p)^7r 


1 — exp 


—)l' 

U+T^ ) \ 


( 21 ) 


where we use +4Z?/T7r — u+ r 2Dj Here the first term corresponds to the pions reaccelerated at the 
shock, and its damping length scale is identical to the distance over which a pion is advected with the fluid 

during its lifetime. The second term represents the pions that are produced from protons in the downstream region 
and simply advected further downward. On the other hand, in the upstream region, Ea. dT^ can be approximated as 


A.- 



1 E>(p)Q,r,o(p) 




exp 




1 D{p)Q^^o{p) 

V. - e 


exp 



( 22 ) 


where we use 



-I- AD/tt^ + u- ^ 2u-. 


B. muon acceleration 

Using the results of the last subsection, we can evaluate the distribution function of muons that are produced from 
the decay of pions. The production spectrum at the source of muons should be given based on the distribution 
function of pions as follows: 


^ _ f-K,±{p/ Vm) _ I f-ir ,±{P / /no-X 

TTvip/^fj.) dp T^ip/^IJ.) 

where R 0.75 is the ratio of the momentum of a muon p to that of a primary pion p^^. Using the method shown 
in the Appendix A, we can solve the muon convection-diffusion equation and derive ffi,±{x,p). The detail of the 
solutions is shown in the Appendix B. 

In the limit of H/uV -C (he., the lifetime of a pion is much longer than the acceleration timescale), the muon 

distribution function at the shock front f^,fi{p) can be approximated as 


f 


7 — a -I- 1 




7 

7 — Q: -|- 1 


1 




2 

TT . 



Qnflip): 

T-jr 


(24) 


where we assume the power-law spectrum of pion production, Qtt,o{p) oc p~°‘. Since D(p) oc p and Ttt oc p, we can 
see that this spectrum Eg. (12411 is proportional to which is similar to the pion distribution function at the shock 

shown in Eq. (12011 . This can be interpreted as follows. As stated in Eg. (12311 . the production spectrum of muons 
is proportional to fn/TTr- Since the lifetime of a pion Tt^ is proportional to p, the production spectrum at the shock 
is proportional to p““+^/p = p““. Injected muons are reaccelerated at the shock and, according 

to the similar discussion to that on pions, the muon spectrum at the shock would become harder than their injected 
spectrum by p^, which comes from the dependence of D{p) on p. We should also note that ^ is proportional to 

(^acc/"Tr)^ accQ7r,0- 
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In a similar way to the pion distribution function, the muon distribution function in the downstream can 

be divided into two components, f^^a.cc{x,p) and ffi^nonaccix^p) as follows: 


ft,,acc{x,p) = //i.oexp - 


^D/t^ - u+ 


2D 


//i,nonacc(^ 5 P) — //i ,0 6 Xp I 


^Ju\+ ^D/t^ - u+ 


2D 


X , 


(25) 


(26) 


and, as in the case of pions, the number of upstream muons, dx^f/^^-(x,p) is subdominant compared to that of 
downstream muons ffj,^+{x,p). 


C. neutrino spectra 


From the pion and muon distribution functions calculated above, the neutrino spectrum can be obtained as follows: 


= I 
= J 

Kip) = I 


fTr{x,p/^„J 

K Mp/U) ’ 

^dnp^ f)j,{x,p/^p^) 

"" K Kp/K) ’ 
Kp/CJ ’ 


(27) 

(28) 
(29) 


where and are the ratios of the energy of a muon neutrino, an anti-muon neutrino, and an electron 

neutrino to that of their primary particles, respectively. Since each lepton produced from the decay of a pion (e, 
and Ve) carries approximately equal energy (i.e., 1/4 of that of the primary pion), we set ~ 0.25, « 0.33 

and ^p^ ~ 0.33. The volume integral should contain the surface area integral on the shocked matter plus the integral 
along the normal direction of the shock. Especially, defining the dynamical timescale tdyn as time for the shock to 
cross the system, the latter integral should be from x ~ — /3_ctdyn to x ~ /3+ctdyn- 

We divide the neutrino energy spectrum into two components according to the decomposition of the pion/muon 
distribution functions shown in Eqs. (Ull), dUl), ([251) and (l26)) : 







,nonacc 



3 47rp^ /7r,acc(2:,p/5pj^ ) 

r.ip/i.,) ’ 

3 47rp^ /tt ,nonacc {X,P/^U^) 

K r^iP/K) 


(30) 

(31) 


and / and , are defined in similar ways. 

i/'^,acc/nonacc i/eiacc/nonacc 

We should also consider neutrino oscillations during the propagation from the source to the Earth. When neutrinos 
propagate over the distances much longer than ^ ficCp/Am^c^ {Am? is the squared mass difference: Ato ^2 — 8.0 x 
10“^ eV^, |Am23| — 2-5 x 10“^ eV^), the observed fluxes of neutrinos $p^ (x = e,p,T) should be described as 


$ 


Prr 


y y i 


(32) 


where U^i is the neutrino mixing matrix and the subscript i represents the mass eigenstate of neutrinos. The matrix 
elements of Uxi can be described by the mixing angles 0i2, 023, and 03i, and the Dirac phase S. Based on [s^, we 
adopt sin^ 0i2 — 0.31, sin^ 023 — 0.39, sin^ 03i ~ 0.024, and S ~ I.Itt. 


III. APPLICATIONS TO LOW-POWER GRBS 

Now we consider long GRBs as neutrino sources. GRBs are thought to produce high-energy neutrinos [l|. In the 
standard model, the emission of long GRBs is believed to be produced by relativistic jets launched when a massive star 
collapses and a stellar-mass black hole is formed. In order for the jet to be observed as a GRB, it should penetrate the 
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stellar envelope, otherwise the jet would stall inside the star and the gamma-ray emission would not be observed [3^ . 
Their prompt emission is often interpreted as synchrotron emission from non-thermal electrons accelerated at internal 
shocks. It is natural to consider the proton acceleration and the associated production of high energy neutrinos via 
PP/Pl interactions [I|. 

However, IceCube gave stringent upper limits on GRBs [13, and has ruled out the typical long GRBs as the 
main source of the observed diffuse neutrino events . 

Instead of ordinary GRBs, we investigate high-energy neutrino production by low-power GRBs such as low- 
luminosity GRBs (LLGRBs) and ultralong GRBs (ULGRBs), which are still not strongly constrained by IceGube. 
In these low-power GRBs the high energy neutrinos may be produced inside the progenitor star . While a jet 

is penetrating in a stellar envelope, it becomes slow and cylindrical by passing through the collimation shock. The 
internal shocks would also occur when there is spatial inhomogeneity in a jet. Murase & loka @ recently investigated 
such high energy neutrino production expected from LLGRBs and ULGRBs [H, ji^l , which have longer durations 
10^ — 10“^ s) and lower luminosity {L^ ^ 10^® — 10^° erg s“^) compared to those of typical long GRBs. It has been 
suggest ed that ultra-long GRBs have bigger progenitors like blue supergiants (BSGs) with radii of ^ 10^^ — 10^^ cm 
[48|-l50|. We apply our model of neutrino production in such GRB jets inside stars, taking into account the secondary- 
acceleration and decay of pions/muons that are produced by shock-accelerated protons via py interactions. The 
internal shocks of GRBs are considered to be mildly-relativistic in the shock rest frame. 

Let us evaluate the important timescales in our model by considering the internal shock scenario of GRBs. When 
two moving shells are ejected with comparable Lorentz factors of order of T from the central engine during the time 
separation At, these shells collide and make an internal shock at the radius r ^ r-i^ ^ T^cAt. Here the magnetic field 
energy density can be estimated as Ub = Lb/ (47rr? cT^), where Lb is the magnetic luminosity. Then we can estimate 
the acceleration timescale face synchrotron cooling timescale U^syn as functions of the energy of a particle, and the 
dynamical timescale fdyn at the shock rest frame as follows: 


tarr — 


D{p) r]ei 


u_ 3ceB/ij_ 


~ 4.4 X 10“^ s ^^biooTeyBiAfnis 


rl/2 02 
^B,47P- 


4^T 


^TT.svn — 


9m^c 


~ 3.0 s 


r?Ae 


2^^ms 


9m;.c 

H' 


Lb,47 £ TV, lOOTeV ’ 
4^7 


~ 0.99 s 


r^Ae 


2^^m.s 


'^B,47£At,100TeV 


^dyn — 


""" / 3 _cr 

= O.IO s r2Af„,,/3z\ 


(33) 


(34) 


(35) 


(36) 


where Si = 100 TeV £4,100 TeV is the energy of a particle * (i = tt or p) at the shock rest frame, r2 = T/IO^, 
Aims = A</(10“^ s), L47 = Lb erg s“\ rriTr ~ 140 MeV and ~ 106 MeV are the masses of a charged pion 
and a muon, respectively. 

From Eqs. dS]) and (jH), a pion can be accelerated at the source before its decay when Gcc < i.e.. 


TyTjAfn 

r 1/2 o 

^B,47P- 


< 4.4 X 10^ 


while a muon can be accelerated before its decay when 

Tyr^Afin, 


rl/2 02 


< 4.7 X 10^. 


(37) 


(38) 


Note that, since both face and Ttt (t^) are proportional to the energy of a pion (a muon), these conditions are 
independent of the energy of particles. Under these conditions, pions (muons) can be accelerated at the shock before 












they decay and therefore their spectra would become harder. On the other hand, we can see that the synchrotron 
cooling timescale would be shorter than the acceleration timescale when the energy Si in the shock rest frame is higher 
than Eifi, where 


£^,o - 2.7 X 10^® 


£^.o - 1-9 X IQi® 


eV 


eV 


rrAtU^p 


(39) 

(40) 


In order to evaluate the timescales of inverse Compton cooling and py interaction, we should give the target photon 
spectrum at the local rest frame. In the case of the internal shock occurring inside a star, the accelerated particles 
mainly interact with photons that are produced in the jet head and escape back from there. Here we estimate the 
spectrum of the target photon field according to the procedure adopted in Q. At the head of the collimated jet, the 
photon temperature Tcj is given as 

feB^cj ~ kB • dcTse) 

« 0.52 keVes;_2iB(V■Y^5(^cj/5)-'/^ (41) 


where L = Lb/eb is the total jet luminosity (es = 0.01eB^_2 is the fraction of the magnetic energy), fcs is the 
Boltzmann constant, ctsb is the Stefan-Boltzmann constant, r^s is the radius where a jet becomes cylindrical through 
the collimation shock, and Fcj is the Lorentz factor of the collimated jet (note that this is different from the Lorentz 
factor of the precollimated jet F). The fraction of photons escaping the collimated jet is /esc ~ (^^cjcrTrcs/Fcj)”^ where 
ricj ~ L/(47rr^gFcjFTOpC^) is the comoving proton number density in the collimated jet, and gt is the Thomson cross 
section. Therefore, the number density of the target photons is given as 


"/escU- 


2rc. 

9.8 X 10^1 cm-3 


(rcj/5) 


- 1/2 


(42) 


where n!() = 167r/(3)(fcBTcj)^/(c/i)^ is the comoving photon number density in the collimated jet, and /(n) is the 
Riemann zeta function. We assume that the escaping photon field has a thermal spectrum, 


dn 

de 




1 


ge/fcBTett _ x’ 


(43) 


with the effective temperature of /cBTeff ~ [(F/2Fcj)/esc]^^^fcB7cj. 

The photomeson production (py interaction) timescale can be evaluated as 


^ 2y2 Xc 


t~ = 

PI 


deaT^{e)^{e)e 


/e/273 


7 -2 

dxx —, 
dx 


(44) 


where y^ = Ep/(iripC^), (T,r(£) is the cross section of pion production as a function of photon energy e in the proton rest 
frame, ^(e) is the average fraction of energy lost from a proton to a pion, and eo = 0.15 GeV is the threshold energy [l|. 
In the following discussion, we use the A resonance approximation: cr,r(£) is approximated to be a function with a peak 
at e = Epeak 0.3 GeV, where cr(£peak) — 5 x 10“^® cm^ with the width of Ae ~ 0.2 GeV, and ^(Epeak) = i-K — 0.2. 

Figures 1 and 2 depict the acceleration timescales, cooling timescales via synchrotron emission and inverse Gompton 
scattering, decay timescales of a pion and a muon, the timescale of py interactions, and the dynamical timescale 
(= rd//3_cF) in the internal shock occurring inside a star expected for an ultralong GRB {L = 10^® erg s“^, eb = 0.01, 
F = 80, At = 10“^ s, [3- = 0.5). We can see that, with the current choice of parameters, the acceleration timescales 
of a pion and a muon are shorter than their lifetimes for arbitrary energy range, and that the decay timescale becomes 
longer than the dynamical timescale above the energy of ^ PeV for pions and ^ 10 TeV for muons (at the shock 
rest frame). Note that, since synchrotron cooling timescale for pions and muons becomes shorter than acceleration 
timescale when the energy of particles is larger than ^ 10 PeV, our formalism is not applicable in the energy range 
above ~ 10 PeV. Note also that the efficiencies of pion/muon production would be suppressed in the energy range 
where the timescale of py interactions is comparable (10^"'^ eV ^ Ei < 10^® eV) or shorter than the acceleration 
timescale. In the current work, this effect is not taken into account. 
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Figure 3 depicts the energy spectra of muon neutrinos and electron neutrinos expected from the internal shock 
inside a progenitor of an ultralong GRB. The flux from reaccelerated pions and muons, Eqs. m and (123, and that 
from advected pions and muons, Eqs. dun and (l26ll . are also shown (with dotted lines and dashed lines, respectively). 
We can see that the electron neutrino flux from advected muons drops above the energy of ^ 100 TeV in the observer 
frame TeV at the shock rest frame). This corresponds to the energy at which the muon decay timescale is equal 
to the dynamical timescale. Above this energy, only part of muons can decay into Ve within the dynamical timescale 
51j . As for the muon neutrino flux from advected particles, it slightly drops around the energy where the electron 
neutrino flux drops because anti-muon neutrinos are generated from the decay of muons fj,'^, and drops again at the 
energy where the pion decay timescale is equal to the dynamical timescale ~ 0.1 PeV for the current parameter 
set) because muon neutrinos are generated from the decay of pions 7r+. We can interpret this behavior as follows. 


From Eqs. (HU), dill), (EZl), (US) and (123), in the limit of ta 
advected particles can be approximated as 


^ Ttt , T'u and tdvn ^ Ttt , Tu, the neutrino fluxes from 


$0 

,nonacc 


(p) + 


,nonacc 

,nonacc 


(p) 

(p) 


V ■ 47r/ [Qn,o(p/^i^J + , 


(45) 

(46) 


while in the limit of face "C tv, T/i and tdyn ^ TV, T/i they can be approximated as 

4’°,,„o„acc(4') - V ■ 4Trp^^-^(tdyn/Tf,)Q^,o(p/^M^^J, (47) 

4>°^.nonaccb) + 4-°^ ,no„acc (p) - ^ ■ dV^dyn [r~^Q^,o(p/^^J + (^^T^ ) "^ (57r,0 (p/Cm?P j] > (48) 

where V is the volume of the merged shell making the internal shock. Here we neglect the contribution from pi- 
ons/muons in the upstream region {fT^/^ _{x,p)) because it is subdominant compared to that from the downstream 
pions/muons. We can easily see that in the latter limit tdyn ^ ''"td'Tij the energy spectra of neutrino fluxes are softer 
than Qtt.o by p^ because the decay timescale Ti is proportional to p. 

On the other hand, the neutrino fluxes from reaccelerated pions/muons increase more as the acceleration timescales 
become longer. Under the condition face "C from Eqs. O, (123, (123, (123 and (123, we can approximate the 

neutrino fluxes from reaccelerated pions/muons as 


4’°e.acc(p) - |l - exp | > 

^°^,acc(p) + 4>°^,acc(p) - Su+Anp^ p{PI11 - CXp | 


+W/m.o(p/CpJ i 1 - exp - 


6,.»'d/r\ 


U+T^ 




(49) 


(50) 


In the high energy limit, where the decay timescales of pions/muons are much longer than the dynamical timescale, 
each of these neutrino fluxes behaves asymptotically as 


‘ Tp. ' ' 

T. 4TTp‘^f^p{p/^^ ) ^ tacc 

$ - U- ^ OC p -, 

T-^ T-jr 




V- 


47rp^fft,o(p/^pJ 


OC 


ba 


' TZ ' fl 


(51) 

(52) 

(53) 


where we use the definition face = D{p)/v?_ and the approximate expressions, Eqs. (123 and (12^ . 

Figure 4 depicts the neutrino flavor ratios as functions of energy expected from the internal shock of ultralong GRBs 
occurring inside progenitors. In addition to the plot for the parameter set used in the previous figures (solid line), we 
show the ratio in the case with longer acceleration timescale for comparison (dashed line). In the usual case, the flavor 
ratio expected from the photomeson process is : 4)°^ =1:2:0, being independent of energy. However, 

when the decay timescale of a muon becomes longer than the dynamical timescale, the flavor ratio is modified because 
the decay timescale of a muon is ^ 100 times longer than that of a pion and only the Ve dux is reduced. On the other 
hand, the acceleration of pions and muons also modifies the flavor ratio, and dominates the neutrino fluxes when the 
acceleration timescale becomes comparable to the dynamical timescale. The flavor ratio becomes constant in the high 
energy limit. We can explain this behavior from Eqs. (123, (123, (E3, (113 and the ratio is determined 

only by the ratio between the acceleration timescale Gcc and the decay timescale of a muon , which is independent 
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of momentum p. More explicitly, when assuming a strong shock (cr = 4, 7 = 4) with Qttp(j>) being proportional to 
p~'^ (i.e., a = 4), we can describe the flavor ratio at the source in the high energy limit as 


*0 1 *0 fa-2ea-l 

(«. (c, 7-«+l + 


*0 — 

I'e 

f-dit + d 

^acc 


~ 0.022^ + 1. (54) 

^acc 

This ratio diverges in the limit of r^/tacc —t 00 , which means that the flavor ratio at the source, 

approaches 0:1:0. Interestingly, we may be able to infer the particle-acceleration timescale from the neutrino flavor 

ratio. 

By using Eq. dal]), we can evaluate the neutrino flavor ratio that would be observed at the Earth, as shown in 
Figure 5. Similar to the flavor ratio at the source, the observed ratio is modified above the energy where the decay 
timescale of a muon becomes longer than the dynamical timescale and is nearly constant in the high energy range. 
The flavor transition occurs over ^ 2 decades in energy. 

We can easily show that, in the limit of r^/tacc —t 00 , the observed flavor ratio high energy 

range converges to ~ 1 : 1.8 : 1.8. This ratio is identical to that shown in Q, in which they investigated the effects of 
synchrotron cooling of pions/muons before their decay on the neutrino flavor ratio. In their work, as in our current 
study, the neutrino flavor ratio at the source is 0 : I : 0 at high energy, but the reason is different. In the model of 
since the lifetime of a muon is longer than that of a pion, muons would suffer from synchrotron cooling more than 
pions. As a result, the flux of electron neutrinos, that are produced from muons, would be suppressed compared to 
the flux of muon neutrinos, that are produced from pions. Therefore, in the high energy range where the synchrotron 
cooling timescale is much shorter than the lifetime of a muon, the flavor ratio at the source can be approximated 
as ~ 0 : 1 : 0. In our study, we show that the flavor ratio would be also modified by the secondary-acceleration 
because pions decay more than muons during the secondary-acceleration, if the acceleration timescale of a pion/muon 
is shorter than their lifetimes. 

On the contrary to the cooling case in which the flavor modification is associated with the spectral softening, the 
neutrino spectra are flat in the high energy range when pions and muons are reaccelerated. This is because the 
secondary-acceleration makes the spectra of primary pions/muons harder by p^ [oc D{p) in Eqs, (I^Dl) and ((Ml) ], while 
the neutrino spectra is softer than the primary spectra by p ^, which is proportional to the decay timescales of primary 
particles. As a result, the neutrino spectra are flat, having the same spectral indices with those of injected primary 
mesons. Therefore, even when the observed flavor ratio of neutrinos in the high energy range converges to 1 : 1.8 : 1.8, 
we can discriminate which process modifies the flavor ratio, cooling or secondary-acceleration, by observing their 
energy spectra. We should note that the flat part of neutrino spectra would have a cutoff at the energy where the 
acceleration timescale is equal to the dynamical timescale because above that energy pions and muons would suffer 
from adiabatic cooling, which is not included in our formulation (see discussion in Sec. 4). 


IV. DISCUSSION AND CONCLUSION 

We investigate the shock acceleration of pions/muons produced by primary protons that are accelerated at the 
shock, and its effects on the observed neutrino flavor ratios. We solve the convection-diffusion equation of pions/muons 
around a shock taking secondary-acceleration and decay into account, and compute the high energy neutrino spectra 
from their decay as well as the energy dependence of the neutrino flavor ratio : 4)^ : We find the following: 

1. When the acceleration timescale is shorter than the decay timescales of a pion and a muon, pions and muons 
are accelerated at the shock before they decay. The resulting distribution function of pions/muons would be divided 
into two components: the component accelerated at the shock and the component advected to the downstream after 
production from protons. The neutrino spectrum of the former component is flat in the high energy range where the 
acceleration timescale becomes comparable to the dynamical timescale of the system. 

2. The flavor ratio of neutrinos at the source, would deviate from 1:2:0, which is expected 

from photomeson interactions, and approaches to 0 : 1 : 0 above the energy at which the decay timescale of a muon 
becomes longer than the dynamical timescale of the shock because only the Ve flux is reduced. The transition width 
of the observed flavor ratio is ^ 2 decades in energy (Eig. 5), which is wider than that in the case of the flavor ratio 
modification W the radiative cooling. Although such a flavor ratio modification by the adiabatic cooling has been 
suggested in [^ , we investigate them using the convection-diffusion equation for the first time. 

3. When the secondary-acceleration is efficient, the neutrino fluxes from shock-reaccelerated pions/muons are 
dominant over the fluxes from non-reaccelerated pions/muons in the high energy range. In this case the flavor ratio 
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would be asymptotically constant (Fig.4). This ratio is determined by the ratio of the lifetime of a muon to its 
acceleration timescale (see Eq. [Ml) . Therefore, from the observed flavor ratio, one can constrain the acceleration 
timescale of cosmic ray particles. 

4. The maximum energy of accelerated particles is determined by the condition face = tdyn, where the energy 
spectra of neutrinos become flat. As a result, the secondary-accelerated component appears as a flat excess above the 
non-reaccelerated component at the highest energy. 

5. When the acceleration timescale is shorter than the lifetime of a muon, the flavor ratio at the source approaches 

to —>■ 0 : 1 : 0 in the high energy range, and the observed flavor ratio approaches to 1 : 1.8 : 1.8. 

This asymptotic ratio is similar to the case where pions/muons are efficiently cooled via synchrotron and/or inverse 
Compton scattering, but the energy spectra of neutrinos are different: the spectra become flat in the high energy 
range when the secondary-acceleration is efficient, while the spectra become soft in the high energy range when the 
synchrotron cooling is efficient. 

6. As for the ratio of f/g to the total v flux, when tacc ^ 1/14 (^ 1/6) in the low energy range and it 

approaches to ^ 0 (^ 1/9) in the case of py {pp) interactions. The ratio of the flux to the total v flux in the high 
energy range can be measured from the j/g interactions at the 6.3 PeV Glashow resonance. 

Our formalism presented in Section II can be applied only when the flow speed at the shock rest frame is non- 
relativistic. In the case when the shock is relativistic, we should use relativistic formulae for shock-acceleration, in 
which the anisotropies in the angular distribution of accelerated particles are taken into account . Recent particle- 
in-cell (PIC) simulations of relativistic shocks have shown that the efficiency of particle acceleration is controlled by 
the magnetization, flow velocity, and field direction [ssl - l^ . and one should take into account these properties when 
discussing secondary-acceleration in relativistic shocks. These issues would be important in future work. 

In the calculations above, we neglect the radiative cooling of pions and muons during the shock acceleration. If 
the energy of pions or muons is higher than in Eqs. dSi) and (HU!) , we should consider the synchrotron cooling in 
deriving the distribution functions of pions and muons. As shown in [3, due to the synchrotron cooling of pions and 
muons, the observed neutrino flavor ratio, $1,^ : is modified from 1 : I : I at low energy to 1 : 1.8 : 1.8 at 

high energy. In this energy range, the energy spectra of neutrinos are softened. This expectation will be confirmed 
by solving the pion/muon transport equations with the energy loss term (e.g. We also neglected the effect of 

matter oscillations (Mekheyev-Smirnov-Wolfenstein effect), which would be important in the case of neutrino emission 
from the GRB jet inside a star because of the high density [I|, These are interesting future works. 


We thank K. Kohri, K. Asano, R. Yamazaki, H. Takami, K. Murase and K. Kashiyama for useful comments. This 
work is supported by the Grants-in-Aid for Scientific Research No. 26287051, 26287051, 24103006, 24000004 and 
26247042 (K.I.). 


Appendix A: General Solutions of the Convection-Diffusion Equation for Decaying Particles 


In this section we discuss how to solve the shock acceleration of charged particles decaying in a finite time such as 
pions and muons. Their convection-diffusion equation is shown in Eq. ©■ The general solution can be described as 


f-lT ni = / + ^ o)> 

I /o°°-I- i7i,+ (x,p) (a: > 0), 


(Al) 


where Gi^±{x'] x,p) are the Green functions of Eq.Q with respect of x, and Hi^±{x,p) are the homogeneous solutions 
of Eq. Q which should be determined by the boundary conditions. 

The Green functions of Eq. Q are given by 


Gi,±(x';x,p) = 


^ul+iD/r, 

1 


exp 


exp 


^ul+4D/ri 

and, under the condition 0, the homogeneous solutions should be 


Ju\Jr4D!Ti^rU^ . ,, 

-2S- (x-x) 


(x > x'), 
(x < x'). 


(A2) 




f -00 dx'Q^ix',p) exp 


\/vF^^-AD7tI—u— , 


2D 


exp 


yj u\_ +4D jri +n_ 

275 • 


(x < 0), 
(A3) 


-C dx'Q^(x\p) exp 3./^ g^p > q). 
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FIG. 1. The acceleration timescale, cooling timescales via synchrotron emission and inverse Compton scattering, and decay 
timescale of charged pions vr'*' in the internal shock occurring inside a progenitor of an ultralong GRB (measured at the shock 
rest frame). The photomeson timescale and dynamical timescale are also shown. Used parameters are Lb = 10“*^ erg s“^, 
r = 80, At = 10"® s, p- = 0.5, and es = 0.01. 


The differential equation for with respect of p is given by the condition (iii) in Eq. (jH]), which can be rewritten 


as 


P~^ = + 19i{p), 


(A4) 


where 7 = 2>aj{a — 1) {a = u_/u+ is the compression ratio), Ai is the numerical factor, being independent of p (since 
both D and Tfi are proportional to p ): 


A; = 


iD 


'1 H-2" ^ I I U 2 ' 2 

TiU_ j \ V O' TiU_ a 


1 AD 1 

+ 


(AS) 


and gi{p) is given by 


gi{p) = — 

U- 


rO I + AD/vi — U- 

dx'Qi{x',p) exp 


One can generally solve Ea. (IA4l) 


as 


2D 


x' \ + J dx'Qi{x' ,p) exp I — 


u\ + AD/n + u+ 
2i4 




x' (.\6) 


(A7) 


Appendix B: Solution of the Convection-Diffusion Equation for Muons 

In this section we describe the solution of the convection-diffusion equation for muons (decaying tertiary particles) 
in detail. The production spectrum of muons at the source (per unit time, per unit spatial volume, and per unit 
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FIG. 2. The same as in Fig. 1, but for muons /r^. 


volume of the momentum space) should be given based on the distribution function of pions as follows: 


^ 1 fTv-ip/ifj.) 


+ 4 :D/Tt^ + U- 




I I — I 

^ I + exp ( -—a; ) , 


D 


^ 1 U,+{p/^t^) + 


Ju\ + 4D/r^ - u+ 




X +a 


/r,b’ 


(Bl) 

(B2) 


where ~ 0.75 is the ratio of the energy of a muon to that of a primary pion, and the functions 9)^a(p) ^-'^el <7^ b(p) 


a^,b'' 


are given as 


'?u,a(p) = — 




f-K,o{p/^f^) - 


_ D{p)Q7v,o{p/^tJ.) 

Cm {D/t^ + (Ctt - C?)M-) 

-D(p)Q^,o(p/Cm) 


Cm {D/t^ + (Ctt - ilW-) 

+ / \ _ /7 i’,o(t/Cm) ^ f /c \ 

Q^,a\P) ~ A Q7r,o(75/^/.i); 

'TV q/x 


^/x,a 


'?M.b(7’) = ^<3 ^.o(t/Cm)- 

S/i 


(B3) 

(B4) 

(B5) 

(B6) 


Note that, since T^(p) is approximately proportional to p, the momentum dependence of Q^(p) is softer than f-Trip)- 
Substituting this Q^{x,p), we can obtain the muon distribution function in the upstream /^^_(x,p) and downstream 
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FIG. 3. The energy flux of (red lines) and t'e (blue lines) expected from a low-power GRB (where the flavor oscillation 

during propagation is not taken into account), normalized to the flux of electron neutrinos at low energy. Used parameters 

are the same as in Fig. 1., and the pion production spectrum is assumed as Q-n,o{jp) oc p~°‘ where a = 4. The muon neutrino flux 
and the electron neutrino flux are divided into two components: those coming from reaccelerated pions and/or muons (dashed 
lines) and those coming from the pions and/or muons advected to the downstream region (dotted lines). In the low energy 
range, the zze flux and flux are dominated by the latter component [oc Qtt.o, see Eg. (1451) and (1461) ]. The Ue, and -\- Up, 

fluxes drop above the energy where the decay timescales of a muon and a pion are equal to the dynamical timescale, being 
proportional to ~ (3,v,ofdyn/v ~ Q,r,ofdyn[l/'r^ +C 2 /tt^], respectively, where C 2 oi oi 10“"' is a constant 

[see Eqs. (1471) and (1481) 1. In the high energy range, the fluxes are dominated by the neutrinos from reaccelerated pions and/or 
muons: the and Up -|- Up fluxes are proportional to ~ Q,r,otLc/TrT> and ~ QTT,o[ti.cc/T-KTp -|- Citacc/T'Tr], respectively, where 
Cl is a constant, the coefficient in front of Tpjt^^^ in Eq. (|54)) [see Eqs. dSU, dsn and (1531) 1. The sum of two components are 
shown by solid lines. Note that, if the energy is higher than 2 x lO'® eV at the observed frame (-^ 3 x lO'® eV at the shock 
rest frame), the spectra would have a cutoff due to the synchrotron cooling of pions and muons (see Figs. 1 and 2), which is 
not taken into account in the current calculation, and therefore the plots above this energy (shown with thin grey lines) would 
be suppressed. Note also that, if the energy is higher than a few times ~ lO'^ eV at the observed frame (a few times 10^® eV 
at the shock rest frame), where the acceleration timescale is longer than the dynamical timescale, the neutrino spectra would 
have a cutoff because there would be no accelerated particles generating neutrinos with such energy. 
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E^, [eV] (shock rest frame) E^, [eV] (shock rest frame) 


FIG. 4. Energy dependence of the flavor ratio of Vfj, + to (left) that of the ratio of to the total neutrino flux 
(right) at the source for low-power GRBs. Used parameter sets {Lb,T^, ^t, (3-) are (10“^^ erg s“^,80,10~® s, 0.5) (solid line) 
and (10"^® erg s“^, 10^,5 x 10“® s, 0.5) (dashed line). As stated in the caption of Fig. 3, in the higher energy range where the 
acceleration timescale is longer than the cooling timescale and/or the dynamical timescale (shown with thin grey lines) the 
ratios would be modified from those shown in these figures. 



Ey [eV] (observed) E^ [eV] (observed) 


FIG. 5. Energy dependence of the flavor ratio of to Vs + Fe (left) and that of the ratio of Vs + to the total neutrino 

flux (right) observed at the Earth for low-power GRBs. Used parameter sets (Lb, T, At, (3-) are (10^^ erg s“^,80,10“® s, 0.5) 
(solid line) and (10^® erg s“^,10^,5 x 10“® s, 0.5) (dashed line). Flavor oscillation during propagation is taken into account. 
As in Figs. 3 and 4, in the high energy range (thin grey lines) the ratios would be modified. 
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ff^,+ ix,p) as 

//J,— ~ I fiJ.,0 


ADq^ 


/j.,a 


Dq 


/j.,b 


(u?. + ID/r^) - - (1 - 


iv?_ + 4Z?/t^ + u_ 

-2:0-^ 




/j,,a 


exp 


+ ^D/Tt^ + u_ 


(m?. +4D/r^) - +4:D/t„ - (1 - \ 


Dq, 


/j.,b 


■ exp 






(B7) 


/ai.+ — I //i.O 


41^0+ 


fj,,a. 


{ul + 4 :D/t^) - ‘iD/r^ + (1 - S,^l)u+Y 


qthT^^ I exp 


\Ju\+ 411/ 


Tfl - u+ 


^Dq^u 


exp - 


\Ju\+ AD/Tt, - u+ 


(m^ +4£)/r^) - {£_^^u\ + AD/t^ + (1 - C^)w+)2 
At the shock front, the muon distribution function should satisfy 


‘iD/i^ 


2D 


X I +9^b^M- 


- -7^/i//x.O + ^ + q^L,hDf_,,h + + q^,hD'^,h) J 


7-0(P) D- 


dp ' ' u 

where A^, and are numerical factors, being independent of p: 


Au = - 

^ 2 


0^,a = 


AD 

1 H-2—I” 4 


1 411 1 


— + ' 2 

O' TnU_ a 


1 + AD/Tfj,u^ ~ 1 + ^Y^^4~~b^O/7vu^ + 1^ 

2 


Dfi^h ~ 


BU = 


dU = 


y/l + AD/t^uI. - (1 - 2^^^^) ' 
2cr 


1 + ADjT^lL^ + 1 + ^Y^l + AD/TTrU^ — 

2cr 


'\ + AD/t^u\ + 1 

One can solve Ea. (IB9l) as 

dp' fp'Y^^ Dip') 


(B8) 

(B9) 

(BIO) 

(Bll) 

(B12) 

(B13) 

(B14) 


U,0ip) ^ (“) (9M.a(p')0^,a + 9;..b(p')0^,b + -l^abOO + a + ^^'.b(p')0+b) ' (^15) 
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